1 Pendahuluan

1.1 Mind Map

1.2 Library

library(tidyverse)
library(lubridate) # date manipulation
library(padr) # TS padding
library(zoo) # imputasi missing value TS
library(fpp) # TS dataset
library(TSstudio) # TS interactive viz

library(forecast) # algoritma forecasting
library(TTR) # SMA function
library(tseries) # adf.test

# supaya semua plot memiliki theme_minimal()
theme_set(theme_minimal())

2 Time Series and Forecasting

Time series data merupakan data yang terurut berdasarkan waktu dan disampel pada interval yang sama.

Forecasting merupakan suatu metode untuk memprediksi/meramalkan data di masa depan

Apa perbedaan Time Series Forecasting dengan Regresi?

Regresi adalah prediksi berdasarkan prediktor atau faktor yang mempengaruhinya. Time series forecasting memprediksi berdasarkan nilai data di masa lalu.

Regression

\[y = \beta_0+\beta_1*x_1+\beta_2*x_2+...+\beta_n*x_n\]

Time Series Forecasting

\[y_t = \beta_0+\beta_1*y_{t-1}+\beta_2*y_{t-2}+...+\beta_n*y_{t-n}\]

Ide utama dalam melakukan forecasting itu adalah korelasi dari data numerik, atau disebut sebagai autokorelasi.

3 Time Series Data

3.1 Characteristics

Time series data: data yang berhubungan dengan waktu dan memiliki interval waktu yang tetap/sama.

Syarat utama data time series:

Data harus urut sesuai periode waktu dari data terlama sampai ke data terbaru
Interval waktunya harus tetap/sama
Tidak boleh data yang terlewat untuk setiap interval

Knowledge Check

Apakah data berikut sudah memenuhi syarat data time series yang baik?

Demand product

df_demand <- data.frame(
  date = ymd(c("2021-5-3", "2021-5-4", "2021-5-6", "2021-5-7")),
  demand = c(29, 79, 41, 88)
  )

df_demand

Kesimpulan: Belum memenuhi karakteristik data TS yang baik, karena ada data yang terlongkap yaitu 5 Mei 2021.

Price product

df_price <- data.frame(
  date = ymd(c("2021-5-16", "2021-5-19", "2021-5-18", "2021-5-17", "2021-5-20", "2021-5-21")),
  price = c(1000, 1001, 1002, 1003, 1004, 1005)
  )

df_price

Kesimpulan: Belum terurut berdasarkan waktu (date)

Perbaikan yang dapat dilakukan sesuai syarat time series:

Mengurutkan data berdasarkan waktu: arrange()

# mengurutkan df_price
df_price %>% 
  arrange(date)

Melakukan padding untuk memastikan interval data sama: pad() dari package padr

Secara default, pad() akan menambal tanggal berdasarkan kolom yang tipe datanya date. Mengisi nilai yang terlewat atau missing (NA), cara yang umum dilakukan dengan package zoo:

na.fill(): mengisi NA dengan sebuah nilai, Gunakan fill="extend" untuk mengisi dengan nilai rata-rata dengan nilai yang missing
na.aggregate(): nilai aggregasi (mean, median)
na.locf(): nilai terakhir sebelum missing

Note: metode untuk mengisi missing value disesuaikan dengan perspektif dari businessnya.

# case: toko tutup di tanggal 5
# kita isi NA dengan nilai 0
df_demand %>% 
  pad() %>% 
  mutate(demand = na.fill(demand, 0))

3.2 Time Series Object

Kita akan menggunakan data emisi CO2 di Indonesia di mana datanya sudah tersimpan dalam folder data_input dengan nama environment_1970f.csv.

# read data
co2 <- read.csv("data_input/environment_1970f.csv")
head(co2)

Data co2 terdiri dari 43 observasi yang mewakili kontribusi gas emisi per tahun terhadap atmosfer Indonesia (43 tahun, 1970-2012). Data ini terdiri dari 7 variabel, yaitu:

year: tahun.
emisi co2 (kt): emisi yang berasal dari pembakaran bahan bakar fosil dan pembuatan semen, termasuk yang dihasilkan selama konsumsi.
emisi co2 (metrik ton per kapita): idem.
emisi metana (kt setara co2): emisi yang berasal dari aktivitas manusia (pertanian) dan dari produksi metana industri.
emisi nitro oksida (ribu metrik ton setara co2): emisi dari pembakaran biomassa pertanian, kegiatan industri, dan pengelolaan ternak.
emisi gas rumah kaca dan lainnya, HFC, PFC dan SF6 (ribu metrik ton setara co2): emisi hasil samping dari hidrofluorokarbon, perfluorokarbon, dan sulfur hexafluoride.
total emisi gas rumah kaca (setara dengan co2): total CO2 tidak termasuk pembakaran biomassa siklus pendek (pembakaran limbah pertanian dan savannah), tetapi termasuk pembakaran biomassa lainnya (kebakaran hutan, pembusukan pasca-pembakaran, kebakaran gambut, dan pembusukan lahan gambut yang dikeringkan), semua sumber CH4 antropogenik, sumber N2O dan gas-F (HFC, PFC, dan SF6).

Dari data co2 ini, kita akan menggunakan 2 kolom:

year untuk menunjukkan waktu
CO2.emissions..metric.tons.per.capita. sebagai nilai yang kita amati untuk membuat object ts

Mengetahui range atau periode waktu data time series, gunakan range()

range(co2$year)

#> [1] 1970 2012

Cek apakah data sudah memenuhi syarat data time series yang baik?

pastikan terurut: ya, sudah terurut
interval waktunya tetap: ya, data tahunan
tidak ada waktu yang terlongkap: ya, semua sudah lengkap

# memastikan apakah tahun sudah terurut dari 1970 sampai 2012
all(co2$year == seq(1970, 2012))

#> [1] TRUE

seq(1970, 2012)

#>  [1] 1970 1971 1972 1973 1974 1975 1976 1977 1978 1979 1980 1981 1982 1983 1984
#> [16] 1985 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999
#> [31] 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012

Mengetahui frequency data yang dimiliki berdasarkan:

data yang disusun per periode apa?
pola yang ingin dilihat apa (harian/mingguan/bulanan/kuartalan/tahunan)?

Berikut adalah contoh cara penentuan frequency:

Jika saya memiliki data harian dan ingin mendapatkan pola mingguan, maka frequency = 7
Jika saya memiliki data harian dan ingin mendapatkan pola bulanan, maka frequency = 7 * 4

Knowledge Check

Tentukan frequency dari data time series berikut:

data total pengunjung cafe per jam, ingin dilihat pola harian, dengan asumsi buka 24 jam maka frequency = 24
data total pengunjung cafe per jam, ingin dilihat pola harian, dengan asumsi hanya buka 12 jam, maka frequency = 12
data total sales product per bulan, ingin dilihat pola tahunan, maka frequency = 12
data total sales product per bulan, ingin dilihat pola kuartalan, maka frequency = 3

Membuat object ts dengan fungsi ts()

Parameter:

data: vector numerik dari data time series
start: periode awal (opsional)
- diisi tahun jika data tahunan
- diisi tahun dan quarter jika data quarterly
- diisi tahun dan bulan jika data bulanan
frequency: pola berulang pada data time series

data tahunan, kemudian kita mau lihat pola tahunan

# membuat object ts
co2_ts <- ts(data = co2$CO2.emissions..metric.tons.per.capita., # vektor numerik
             start = 1970, # tahun awal
             frequency = 1) 
co2_ts

#> Time Series:
#> Start = 1970 
#> End = 2012 
#> Frequency = 1 
#>  [1] 0.3119519 0.3306215 0.3580080 0.3954703 0.4021567 0.4128050 0.4612390
#>  [8] 0.6002978 0.6677802 0.6601457 0.6426495 0.6634071 0.6822242 0.6640976
#> [15] 0.6944021 0.7347680 0.7229173 0.7184145 0.7552095 0.7348064 0.8243417
#> [22] 0.9735380 1.0788747 1.1452296 1.1416286 1.1420774 1.2669930 1.3738789
#> [29] 1.0218517 1.1599925 1.2452416 1.3748183 1.4102338 1.4364045 1.5098982
#> [36] 1.5084806 1.5015768 1.6118554 1.7638951 1.8651654 1.7679079 2.3335854
#> [43] 2.4089201

Visualisasi object ts

Gunakan autoplot() dari package forecast, yang mengembalikan object ggplot:

library(forecast)
co2_ts %>% 
  autoplot()

Misalkan ingin cek pola pada tahun 1990 - 2000 saja, gunakan window(start, end)

window: melakukan subsetting pada object time series

co2_ts %>% 
  window(start = 1990, end = 2000) %>% 
  autoplot()

Insight: Pergerakan emisi gas co2 di Indonesia secara umum naik, tapi ada yang menarik di tahun 1997-1998. Terjadi peningkatan emisi gas yang amat tinggi pada tahun 1997 yang mungkin disebabkan oleh Kebakaran Hutan 1997

Diskusi: Objek co2_ts berasal dari data yang direkam tahunan dan kita mengatur pola tahunan (frequency = 1). Apakah kita bisa menganalisis pola berulang (musimannya)?

Jawaban: tidak bisa menganalisis pola berulang apabila frequency = 1

Note: pemilihan frequency umumnya menggunakan ukuran waktu 1 level di atasnya, atau yang lebih atas lagi. (data harian dengan seasonality mingguan / bulanan / dst.)

3.3 Exploring `nybirth` Data

Read data nybirth.csv

birth <- read.csv("data_input/nybirth.csv")
glimpse(birth)

#> Rows: 168
#> Columns: 2
#> $ date   <chr> "1946-01-01", "1946-02-01", "1946-03-01", "1946-04-01", "1946-0~
#> $ births <dbl> 26.663, 23.598, 26.931, 24.740, 25.806, 24.364, 24.477, 23.901,~

Data di atas merupakan data persentase kelahiran di New York per bulan. Terdiri dari:

date: tanggal saat dilakukan pencatatan persentase kelahiran
births: persentase kelahiran

Cek class data birth:

class(birth)

#> [1] "data.frame"

# masih data frame

Lakukan data preprocessing:

menyesuaikan tipe data date
ekstrak informasi bulan dari date untuk keperluan visualisasi nantinya

# gunakan fungsi dari package lubridate
birth_clean <- birth %>% 
  mutate(date = ymd(date),
         month = month(date, label = T))
birth_clean

Mengetahui range atau periode waktu data birth

range(birth_clean$date)

#> [1] "1946-01-01" "1959-12-01"

Cek apakah data sudah memenuhi syarat data time series yang baik?

data harus urut
interval tetap (tidak ada yang terlongkap)
tidak ada missing value

birth_clean %>% 
  arrange(date) %>% 
  pad() %>% 
  anyNA()

#> [1] FALSE

Kesimpulan: data birth_clean sudah memenuhi ketiga karakteristik data TS, karena setelah diurutkan dan dilakukan padding, tidak ada missing value

Eksplorasi data

Visualisasi data sebelum mengubah menjadi object ts untuk mengetahui pola berulang pada data birth

ggplot(data = birth_clean, aes(x = date, y = births)) +
  geom_point() +
  geom_line()

ggplot(data = birth_clean, aes(x = date, y = births)) +
  geom_point() +
  geom_point(data = birth_clean %>% 
               filter(month %in% c("Jan", "Feb", "Jul")),
             aes(col = month)) +
  geom_line() +
  scale_color_manual(values = c("red", "blue", "yellow"))

Insight: Bagaimana pola dari data birth? ada pola yang berulang setiap 12 bulan (pola tahunan) karena setiap Feb tingkat kelahiran cenderung rendah, sedangkan tiap Jul cenderung tinggi.

Membuat object time series

birth_ts <- ts(data = birth_clean$births,
               frequency = 12, # data bulanan, pola tahunan
               start = c(1946, 1)) # opsional, data awal = 1946 Jan
birth_ts

#>         Jan    Feb    Mar    Apr    May    Jun    Jul    Aug    Sep    Oct
#> 1946 26.663 23.598 26.931 24.740 25.806 24.364 24.477 23.901 23.175 23.227
#> 1947 21.439 21.089 23.709 21.669 21.752 20.761 23.479 23.824 23.105 23.110
#> 1948 21.937 20.035 23.590 21.672 22.222 22.123 23.950 23.504 22.238 23.142
#> 1949 21.548 20.000 22.424 20.615 21.761 22.874 24.104 23.748 23.262 22.907
#> 1950 22.604 20.894 24.677 23.673 25.320 23.583 24.671 24.454 24.122 24.252
#> 1951 23.287 23.049 25.076 24.037 24.430 24.667 26.451 25.618 25.014 25.110
#> 1952 23.798 22.270 24.775 22.646 23.988 24.737 26.276 25.816 25.210 25.199
#> 1953 24.364 22.644 25.565 24.062 25.431 24.635 27.009 26.606 26.268 26.462
#> 1954 24.657 23.304 26.982 26.199 27.210 26.122 26.706 26.878 26.152 26.379
#> 1955 24.990 24.239 26.721 23.475 24.767 26.219 28.361 28.599 27.914 27.784
#> 1956 26.217 24.218 27.914 26.975 28.527 27.139 28.982 28.169 28.056 29.136
#> 1957 26.589 24.848 27.543 26.896 28.878 27.390 28.065 28.141 29.048 28.484
#> 1958 27.132 24.924 28.963 26.589 27.931 28.009 29.229 28.759 28.405 27.945
#> 1959 26.076 25.286 27.660 25.951 26.398 25.565 28.865 30.000 29.261 29.012
#>         Nov    Dec
#> 1946 21.672 21.870
#> 1947 21.759 22.073
#> 1948 21.059 21.573
#> 1949 21.519 22.025
#> 1950 22.084 22.991
#> 1951 22.964 23.981
#> 1952 23.162 24.707
#> 1953 25.246 25.180
#> 1954 24.712 25.688
#> 1955 25.693 26.881
#> 1956 26.291 26.987
#> 1957 26.634 27.735
#> 1958 25.912 26.619
#> 1959 26.992 27.897

Visualisasi object ts:

Misal kita tertarik untuk visualisasi birth dari Jan 1950 sampai Des 1955 dengan autoplot()

birth_ts %>% 
  window(start = c(1950, 1), end = c(1955, 12)) %>% 
  autoplot()

Opsional: Gunakan ts_plot() dari package TSstudio untuk visualisasi interaktif:

birth_ts %>% 
  ts_plot()

3.4 Decomposition

Decomposition adalah suatu tahapan dalam analisis time series untuk menguraikan data menjadi beberapa komponen dalam time series data, yaitu:

Trend: pola data secara umum, kecenderungan untuk naik atau turun
Seasonal: pola musiman yang membentuk pola berulang pada periode waktu yang tetap
Error/Remainder/Random: pola yang tidak dapat ditangkap dalam trend dan seasonal

Untuk dapat menguraikan object time series kita menjadi 3 komponen tersebut, kita dapat menggunakan fungsi decompose().

birth_dc <- birth_ts %>% 
  decompose()

Visualisasi hasil decompose:

birth_dc %>% 
  autoplot()

Notes: Jika pada hasil decompose, trend masih membentuk sebuah pola maka dapat dicurigai masih ada seasonality yang belum ditangkap. Seharusnya trend cenderung naik atau cenderung turun secara smooth. Penyebabnya:

frequency yang kita tetapkan belum tepat, atau
terdapat multiseasonality pada data (keep untuk materi additional)

# coba buat ulang object ts dari birth_clean tapi frequency = 4
ts(data = birth_clean$births, frequency = 4) %>% 
  decompose() %>% 
  autoplot()

Notes: Object time series dengan frequency = 1, tidak bisa dibuat decomposenya

# co2_ts %>% 
#  decompose()

# time series has no or less than 2 periods

3.5 Additive and Multiplicative

Terdapat 2 jenis model pada data time series, yaitu:

Model Additive: Model time series yang memiliki varians konstan mengikuti trend dan seasonalnya

\[Y_t = T_t + S_t + E_t\]

Data = Trend + Seasonal + Error

Contoh time series additive:

birth_ts %>%
  autoplot()

Model Multiplicative: Model time series yang memiliki varians semakin tinggi/rendah mengikuti trend dan seasonal yang ada

\[Y_t = T_t * S_t * E_t\]

Data = Trend * Seasonal * Error

Contoh time series multiplicative:

AirPassengers %>% 
  autoplot()

3.5.1 (Optional) Perhitungan manual decompose Additive

Jika kita perhatikan lagi pada data birth_ts memiliki pola additive, karena varians dari polanya tetap atau konstan. Secara default, fungsi decompose() memiliki type = "additive".

Melakukan inspect komponen time series pada data birth_dc

1. Trend

Trend diperoleh dari hasil perhitungan center moving average (CMA). Tujuan utamanya untuk smoothing data sehingga diperoleh trend yang cenderung naik atau turun.

Berikut pendekatan manual untuk perhitungan trend, menggunakan fungsi ma(). Parameter:

x: object TS
order: berapa banyak data yang dilibatkan, menyesuaikan frequency
centre = T: menggunakan center moving average

# perhitungan manual, di mana order = frequency pada object TS
ma(birth_ts, order = 12, centre = T) %>% head(24)

#>           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
#> 1946       NA       NA       NA       NA       NA       NA 23.98433 23.66213
#> 1947 22.35350 22.30871 22.30258 22.29479 22.29354 22.30562 22.33483 22.31167
#>           Sep      Oct      Nov      Dec
#> 1946 23.42333 23.16112 22.86425 22.54521
#> 1947 22.26279 22.25796 22.27767 22.35400

Bandingkan dengan $trend dari birth_dc:

# trend dari hasil decompose()
birth_dc$trend %>% head(24)

#>           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
#> 1946       NA       NA       NA       NA       NA       NA 23.98433 23.66213
#> 1947 22.35350 22.30871 22.30258 22.29479 22.29354 22.30562 22.33483 22.31167
#>           Sep      Oct      Nov      Dec
#> 1946 23.42333 23.16112 22.86425 22.54521
#> 1947 22.26279 22.25796 22.27767 22.35400

# hasilnya sama persis dengan perhitungan manual

2. Seasonal

Additive: Data = Trend + Seasonal + Error Data de-trend: Data - Trend = Seasonal + Error

Data - Trend = Seasonal + Error

birth_detrend <- birth_ts - birth_dc$trend

Pendekatan manual

# mean tiap bulan
mean_month <- birth_detrend %>%
  matrix(ncol = 12, byrow = T) %>% 
  colMeans(na.rm = T)
  
# mean global
mean_global <- mean(mean_month)

# nilai seasonality
mean_month - mean_global

#>  [1] -0.6771947 -2.0829607  0.8625232 -0.8016787  0.2516514 -0.1532556
#>  [7]  1.4560457  1.1645938  0.6916162  0.7752444 -1.1097652 -0.3768197

Bandingkan dengan $seasonal dari birth_dc:

# seasonal dari hasil decompose()
birth_dc$seasonal %>% head(12)

#>             Jan        Feb        Mar        Apr        May        Jun
#> 1946 -0.6771947 -2.0829607  0.8625232 -0.8016787  0.2516514 -0.1532556
#>             Jul        Aug        Sep        Oct        Nov        Dec
#> 1946  1.4560457  1.1645938  0.6916162  0.7752444 -1.1097652 -0.3768197

# hasilnya sama persis dengan perhitungan manual

3. Error

Additive: Data = Trend + Seasonal + Error Additive error: Error = Data - Trend - Seasonal

(birth_ts - birth_dc$trend - birth_dc$seasonal) %>% 
  autoplot()

Bandingkan dengan $random dari birth_dc:

# error dari hasil decompose()
birth_dc$random %>% 
  autoplot()

# hasilnya sama persis dengan perhitungan manual

Inti dari perhitungan manual:

Trend: menggunakan Center Moving Average
Seasonality: Mean berdasarkan masing-masing periode (bulanan)
Error: nilai sisanya, apabila additive yaitu data - trend - seasonality

Istilah:

Remainder: dalam visualisasi
Random: elemen dalam list $random
Error: dalam pemodelan forecasting

3.5.2 Multiplicative Time Series

Jika kita memiliki pola multiplicative pada data time series dan ingin membuat decomposenya cukup menambahkan parameter type = "multiplicative" pada fungsi decompose().

Ketika kita menemukan pola data kita mengandung multiplicative:

Cara 1: transformasi data multiplicative menjadi additive dengan fungsi log. Setelah memperoleh hasil forecast kita dapat mengembalikan nilainya dengan exp.

AirPassengers %>% 
  autoplot()

# transformasi log
AirPassengers %>% 
  log() %>% 
  autoplot()

(Opsional) Sifat logaritma: perkalian menjadi penjumlahan

\[y = T * S * E\] -> multiplicative \[log(y) = log(T * S * E)\] \[log(y) = log(T) + log(S) + log(E)\] -> additive

Cara 2: Tetap menggunakan model multiplicative, kemudian nanti hasil dibandingkan dengan memilih model dengan error yang paling kecil.

Parameter type dalam fungsi decompose(), secara default type = "additive"

air_decom <- AirPassengers %>% 
  decompose(type = "multiplicative")

air_decom %>% 
  autoplot()

3.5.3 (Opsional) Perhitungan manual decompose Multiplicative

Trend

Menggunakan Center Moving Average (CMA)

ma(x = air_decom$x, order=12, centre = T) %>% head(12)

#>           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
#> 1949       NA       NA       NA       NA       NA       NA 126.7917 127.2500
#>           Sep      Oct      Nov      Dec
#> 1949 127.9583 128.5833 129.0000 129.7500

Bandingkan dengan trend di air_decom:

air_decom$trend %>% head(12)

#>           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
#> 1949       NA       NA       NA       NA       NA       NA 126.7917 127.2500
#>           Sep      Oct      Nov      Dec
#> 1949 127.9583 128.5833 129.0000 129.7500

# hasilnya sama persis dengan perhitungan manual

Seasonality

Multiplicative: Data = Trend * Seasonal * Error De-trend: Data / Trend = Seasonal * Error

seasxerr <- air_decom$x /air_decom$trend

# mean of each month
mean_month <- seasxerr %>%
  matrix(ncol = 12, byrow = T) %>% 
  colMeans(na.rm = T)

# mean global
mean_global <- mean(mean_month)

# Seasonality Calculation
mean_month / mean_global

#>  [1] 0.9102304 0.8836253 1.0073663 0.9759060 0.9813780 1.1127758 1.2265555
#>  [8] 1.2199110 1.0604919 0.9217572 0.8011781 0.8988244

Bandingkan dengan seasonal di air_decom:

air_decom$seasonal %>% head(12)

#>            Jan       Feb       Mar       Apr       May       Jun       Jul
#> 1949 0.9102304 0.8836253 1.0073663 0.9759060 0.9813780 1.1127758 1.2265555
#>            Aug       Sep       Oct       Nov       Dec
#> 1949 1.2199110 1.0604919 0.9217572 0.8011781 0.8988244

# hasilnya sama persis dengan perhitungan manual

Error

Multiplicative: Data = Trend * Seasonal * Error Multiplicative Error: Error = Data / Trend / Seasonal

(air_decom$x/air_decom$trend/air_decom$seasonal) %>% 
  autoplot()

Bandingkan dengan random di air_decom:

air_decom$random %>% 
  autoplot()

# hasilnya sama persis dengan perhitungan manual

3.6 Seasonality Analysis

Seasonality analysis membantu kita mengetahui di waktu mana saja yang nilai datanya tinggi/rendah pada periode seasonal yang kita amati. Kita bisa menggunakan informasi $seasonal untuk membuat visualisasi seasonal bar plot:

# analisis seasonal dengan bar plot

birth_clean %>% 
  mutate(seasonal = birth_dc$seasonal) %>% # mengekstrak nilai seasonal
  distinct(month, seasonal) %>% # mengambil baris yang unique saja
  ggplot(aes(x = month, y = seasonal)) + # melakukan visualisasi
  geom_col()

3.7 Seasonality Adjusted

Seasonality adjusted adalah data time series yang sudah dibuang efek seasonal nya. Umumnya digunakan untuk lebih mudah mendeteksi error/kejadian luar biasa/anomali dari data (tidak terganggu efek seasonal). Hal ini untuk kebutuhan exploratory data saja.

Berikut contoh data birth_ts yang sudah dibuang efek seasonalnya:

Model additive:

Data = Trend + Seasonality + Error Data - Seasonality = Trend + Error

# $x = data asli
# $seasonal = seasonality
(birth_dc$x - birth_dc$seasonal) %>%
  autoplot() %>% 
  plotly::ggplotly()

Visualisasi berikut menunjukkan line plot untuk Trend sekaligus Seasonal Adjusted data. Dengan seperti ini, kita dapat dengan mudah mendeteksi anomali (titik-titik data yang jauh dari trend nya).

Credit to Pak Bayu Andrianto :)

birth_dc$trend %>% 
  autoplot(series = "Trend") +
  autolayer(birth_dc$x - birth_dc$seasonal, series = "Seasonal adjusted")

4 📝 SUMMARY DAY 1

Perbedaan regresi dengan forecasting
- regresi: predictornya menggunakan variable lainnya (independent)
- forcasting: predictornya menggunakan variable target itu sendiri di waktu lalu
Tahap pengolahan data time series:

Memastikan data sudah memenuhi syarat:

Tidak ada missing value (NA) -> imputasi
Data urut berdasarkan tanggalnya -> sorting arrange()
Interval waktu tetap, tidak boleh ada yang terlongkap -> padding pad()

Membuat ts object dengan fungsi ts(). Parameter:

data: vektor numerik
start: tanggal awal
frequency: pola berulang yang ingin ditangkap

Visualisasi object ts dengan autoplot():

Mengetahui apakah additive atau multiplicative
Additive: variansi konstan/tetap
Multiplicative: variansi berubah (besar/kecil)

Melakukan decompose dengan decompose() menjadi 3 komponen:

trend: naik/turun
seasonality: pola musiman/berulang
error: pola yang belum tertangkap oleh trend dan seasonality

Tujuan melakukan decomposition: - insight dari trend dan seasonal - mengetahui apakah frequency yang kita set sudah tepat atau belum.

DAY 2

5 Forecasting

Metode:

Simple Moving Average (SMA)
Exponential Smoothing (ets = error trend seasonal)
- Simple/Single Exponential Smoothing (SES)
- Double Exponential Smoothing (Holt Exponential Smoothing)
- Triple Exponential Smoothing (Holt Winters Exponential Smoothing)
Auto Regressive Integrated Moving Average (ARIMA)
Seasonal Auto Regressive Integrated Moving Average (SARIMA)
Additional: Seasonal Trend with Loess Model (STLM)

5.1 Simple Moving Average (SMA)

Metode ini menggunakan rataan bergerak (moving average) untuk melakukan smoothing data. Karena menggunakan rataan, maka bobot yang digunakan sama untuk setiap observasi di masa lalu.

\[y_t = \frac{y_{t-1}+y_{t-2}+...+y_{t-n}}{n}\]

SMA tepat digunakan untuk data yang tidak mengandung trend dan tidak ada seasonal

trend: pola data secara umum apakah naik/turun seasonal: pola berulang

y <- c(3, 4, 5, 6)
mean(y)

#> [1] 4.5

(3+4+5+6)/4

#> [1] 4.5

1/4 * 3 + 1/4 * 4 + 1/4 * 5 + 1/4 * 6

#> [1] 4.5

SMA: berdasarkan rata-rata, dan tiap datanya memiliki bobot yang sama

Import data

Kita menggunakan data curah hujan tahunan dari tahun 1813-1912 (100 tahun) di London, England.

library(TTR)
rain <- scan("data_input/precip1.dat", skip = 1)
head(rain)

#> [1] 23.56 26.07 21.86 31.24 23.65 23.88

Note: skip = 1 untuk melongkap 1 observasi pada file .dat, karena ada tulisan “Total”.

Buat objek time series

data tahunan -> pola tahunan

rain_ts <- ts(data = rain,
              start = 1813, # tahun awal
              frequency = 1)
rain_ts

#> Time Series:
#> Start = 1813 
#> End = 1912 
#> Frequency = 1 
#>   [1] 23.56 26.07 21.86 31.24 23.65 23.88 26.41 22.67 31.69 23.86 24.11 32.43
#>  [13] 23.26 22.57 23.00 27.88 25.32 25.08 27.76 19.82 24.78 20.12 24.34 27.42
#>  [25] 19.44 21.63 27.49 19.43 31.13 23.09 25.85 22.65 22.75 26.36 17.70 29.81
#>  [37] 22.93 19.22 20.63 35.34 25.89 18.65 23.06 22.21 22.18 18.77 28.21 32.24
#>  [49] 22.27 27.57 21.59 16.93 29.48 31.60 26.25 23.40 25.42 21.32 25.02 33.86
#>  [61] 22.67 18.82 28.44 26.16 28.17 34.08 33.82 30.28 27.92 27.14 24.40 20.35
#>  [73] 26.64 27.01 19.21 27.74 23.85 21.23 28.15 22.61 19.80 27.94 21.47 23.52
#>  [85] 22.86 17.69 22.54 23.28 22.17 20.84 38.10 20.65 22.97 24.26 23.01 23.67
#>  [97] 26.75 25.36 24.79 27.88

Apakah data rain dapat di-decompose? Karena frequency = 1, tidak bisa dilakukan decompose

Visualisasikan data

rain_ts %>% 
  autoplot()

Istilah: Data curah hujan tidak memiliki trend dan seasonal, atau bisa kita sebut sebagai data stasioner.

Fitting Model SMA dengan ordo 3

Fungsi yang digunakan untuk fitting model dengan metode SMA adalah SMA() dari package TTR dengan parameter:

x: objek time series
n: ordo

rain_sma3 <- SMA(rain_ts, n = 3) # melihat 3 data ke belakang
rain_sma3

#> Time Series:
#> Start = 1813 
#> End = 1912 
#> Frequency = 1 
#>   [1]       NA       NA 23.83000 26.39000 25.58333 26.25667 24.64667 24.32000
#>   [9] 26.92333 26.07333 26.55333 26.80000 26.60000 26.08667 22.94333 24.48333
#>  [17] 25.40000 26.09333 26.05333 24.22000 24.12000 21.57333 23.08000 23.96000
#>  [25] 23.73333 22.83000 22.85333 22.85000 26.01667 24.55000 26.69000 23.86333
#>  [33] 23.75000 23.92000 22.27000 24.62333 23.48000 23.98667 20.92667 25.06333
#>  [41] 27.28667 26.62667 22.53333 21.30667 22.48333 21.05333 23.05333 26.40667
#>  [49] 27.57333 27.36000 23.81000 22.03000 22.66667 26.00333 29.11000 27.08333
#>  [57] 25.02333 23.38000 23.92000 26.73333 27.18333 25.11667 23.31000 24.47333
#>  [65] 27.59000 29.47000 32.02333 32.72667 30.67333 28.44667 26.48667 23.96333
#>  [73] 23.79667 24.66667 24.28667 24.65333 23.60000 24.27333 24.41000 23.99667
#>  [81] 23.52000 23.45000 23.07000 24.31000 22.61667 21.35667 21.03000 21.17000
#>  [89] 22.66333 22.09667 27.03667 26.53000 27.24000 22.62667 23.41333 23.64667
#>  [97] 24.47667 25.26000 25.63333 26.01000

Visualisasi data actual dan hasil smoothing SMA dengan autolayer()

rain_ts %>% 
  autoplot(series = "Rain") +
  autolayer(rain_sma3, series = "SMA(3)")

# series = label di legend

Knowledge Check: Misal kita gunakan SMA3 untuk forecast tahun 1913, maka nilainya adalah 26.01

# manual
tail(rain_ts, 3) %>% mean()

#> [1] 26.01

# prediksi nilai di masa depan
library(forecast)
forecast(rain_sma3, h = 1)

#>      Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
#> 1913       26.00996 23.68437 28.33555 22.45328 29.56665

Kelemahan SMA:

Hanya n data ke belakang
Pembobotan yang sama pada tiap waktunya, baik untuk nilai di waktu terbaru maupun di masa lalu

5.2 Exponential Smoothing

Pada exponential smoothing, bobot yang digunakan menurun secara eksponensial. Data yang lebih baru akan diberikan bobot yang lebih besar dibandingkan data yang sudah lama. Metode exponential smoothing terbagi menjadi 3 yaitu:

Single/Simple Exponential Smooting/SES: no trend, no seasonal
Double Exponential Smooting/Holt Exponential Smoothing: ada trend, no seasonal
Triple Exponential Smooting/Holt-Winters Exponential Smoothing: ada trend, ada seasonal

Fungsi yang digunakan untuk membuat model exponential smoothing, yaitu:

HoltWinters(), parameter yang digunakan, yaitu:

x: objek time series
secara default parameter alpha, beta, dan gamma adalah NULL. Parameter harus diubah menjadi FALSE apabila komponen tersebut tidak ada.
- alpha = smoothing untuk error, nilainya antara 0 - 1
- beta = smoothing untuk trend, nilainya antara 0 - 1
- gamma = smoothing untuk seasonal, nilainya antara 0 - 1
- seasonality = jenis model, apakah additive/multiplicative

ets -> bobot abg

SES: HoltWinters(x, beta = F, gamma = F)

\[y_{t+1} = \alpha y_{t} + \alpha (1-\alpha)y_{t-1}+\alpha(1-\alpha)^2y_{t-2}+...\]

alpha mendekati 1 artinya observasi yang paling baru diberikan bobot jauh lebih tinggi dibandingkan data di waktu yang lama
alpha mendekati 0 artinya observasi yang paling baru tetap memiliki bobot yang lebih tinggi namun tidak berbeda jauh dibandingkan dengan bobot di waktu yang lama

Double/Holt Exponential Smoothing: HoltWinters(x, gamma = F)

\[y_{t+h|t} = \ell_t +hb_t\] \[\ell_t = \alpha y_t+(1-\alpha)(\ell_{t-1}+b_{t-1})\] \[b_t = \beta^{*}(\ell-\ell_{t-1}) +(1-\beta^{*})b_{t-1}\]

Triple/Holt Winters Exponential Smoothing: HoltWinters(x)

\[y_{t+h|t} = \ell_t +hb_t+s_{t+h-m(k+1)}\] \[\ell_t = \alpha y_t+(1-\alpha)(\ell_{t-1}+b_{t-1})\] \[b_t = \beta^{*}(\ell-\ell_{t-1}) +(1-\beta^{*})b_{t-1}\] \[s_t = \gamma(y_t-\ell_{t-1}-b_{t-1})+(1-\gamma)s_{t-m}\]

5.2.1 Simple Exponential Smoothing (SES)

SES baik digunakan untuk data time series yang tidak mengandung trend dan tidak seasonal. SES akan melakukan smoothing terhadap komponen error/random.

Kita akan menggunakan data curah hujan untuk melakukan fitting model dan forecasting menggunakan metode SES.

rain_ts %>% 
  autoplot()

Cross validation: train-test splitting

splitting tidak boleh random, harus terurut
data test harus data terbaru

# 20 tahun terakhir untuk test
rain_test <- tail(rain_ts, 20)

# 80 tahun untuk train (sisanya)
rain_train <- head(rain_ts, -20) # ekuivalen dengan head(rain_ts, 80)

Fitting Model dengan fungsi HoltWinters()

ets => parameter abg

no trend, no seasonal -> beta = F, gamma = F

rain_ses <- HoltWinters(rain_train, beta = F, gamma = F)
rain_ses

#> Holt-Winters exponential smoothing without trend and without seasonal component.
#> 
#> Call:
#> HoltWinters(x = rain_train, beta = F, gamma = F)
#> 
#> Smoothing parameters:
#>  alpha: 0.03595009
#>  beta : FALSE
#>  gamma: FALSE
#> 
#> Coefficients:
#>       [,1]
#> a 25.24619

Note: secara default fungsi HoltWinters() akan mencari parameter smoothing yang menurutnya optimal -> meminimalisir error di data train. “Optimal” itu belum tentu menghasilkan the best performance untuk di data test.

Apabila ingin set nilai smoothing sesuai keinginan, misal alpha di nilai 0.8 maka: HoltWinters(rain_train, alpha = 0.8, beta = F, gamma = F)

Forecasting dengan fungsi forecast()

rain_forecast <- forecast(rain_ses, h = 20) # panjang data test

rain_forecast %>% 
  autoplot()

Parameter: h = horizon (berapa waktu kedepan yang ingin di-forecast?)

Output:

Point Forecast ($mean): nilai forecasting di masa depan
Lo 80: batas bawah interval forecasting dengan confidence level 80%
Hi 80: batas atas interval forecasting dengan confidence level 80%
Lo 95: batas bawah interval forecasting dengan confidence level 95%
Hi 95: batas atas interval forecasting dengan confidence level 95%

Visualisasi data actual dengan hasil forecasting SES

rain_ts %>% 
  autoplot(series = "Actual") +
  autolayer(rain_forecast$fitted, series = "Train") +
  autolayer(rain_forecast$mean, series = "Test")

Note:

$fitted untuk hasil forecast di data train
$mean untuk hasil forecast di data test

Evaluasi model dengan fungsi accuracy() dari library forecast

Parameter: object_forecast dan object_ts_test

accuracy(rain_forecast, rain_test)

#>                      ME     RMSE      MAE        MPE    MAPE      MASE
#> Training set  0.5937162 4.310376 3.431055 -0.4116781 13.7016 0.6794168
#> Test set     -1.2686879 4.280286 3.248594 -7.8305099 13.7201 0.6432859
#>                     ACF1 Theil's U
#> Training set -0.05786692        NA
#> Test set     -0.24464714 0.7232343

Train: 13.7% Test: 13.72%

Metric evaluasi:

ME: Mean Error
RMSE: Root Mean Squared Error
MAE: Mean Absolute Error
MPE: Mean Percentage Error
MAPE: Mean Absolute Percentage Error
MASE: Mean Absolute Scaled Error
ACF1: Autocorrelation of errors at lag 1.
Theil’s U: RMSE divided by RMSE of naive model $y_t = y_{t-1}$. Idealnya di bawah 1

5.2.2 Double Exponential Smoothing / Holt Exponential Smoothing

Holt digunakan untuk data time series yang mengandung trend, namun tidak seasonal. Holt akan melakukan smoothing terhadap komponen error/random dan trend.

Kita akan menggunakan data kontribusi emisi gas CO2 di Indonesia untuk melakukan fitting model dan forecasting menggunakan metode holt.

co2_ts %>% 
  autoplot()

Note: ingat, karena frequency = 1 kita tidak dapat melakukan decompose()

Cross validation

# 3 tahun terakhir untuk test
co2_test <- tail(co2_ts, 3)

# 40 tahun untuk train (sisanya)
co2_train <- head(co2_ts, -3)

Fitting Model Double Exponential Smoothing (DES)

ets => parameter abg

ada trend, tidak ada seasonal -> gamma = F

co2_des <- HoltWinters(co2_train, gamma = F)
co2_des

#> Holt-Winters exponential smoothing with trend and without seasonal component.
#> 
#> Call:
#> HoltWinters(x = co2_train, gamma = F)
#> 
#> Smoothing parameters:
#>  alpha: 0.9035708
#>  beta : 0.03548442
#>  gamma: FALSE
#> 
#> Coefficients:
#>         [,1]
#> a 1.85771545
#> b 0.03875735

Hati-hati: Jangan set nilai parameter smoothing menjadi TRUE. Contoh kalau beta = T artinya beta = 1

Forecasting

co2_forecast <- forecast(co2_des, h = 3) # panjang data test

Visualisasi data actual dengan hasil forecasting

co2_ts %>% 
  autoplot(series = "Actual") +
  autolayer(co2_forecast$fitted, series = "Train") +
  autolayer(co2_forecast$mean, series = "Test")

Evaluasi model

forecast::accuracy(co2_forecast, co2_test)

#>                      ME       RMSE        MAE      MPE      MAPE      MASE
#> Training set 0.01648723 0.08542822 0.05645794 1.323059  5.813989 0.9019799
#> Test set     0.23490765 0.34851219 0.32061760 9.284487 14.132587 5.1222313
#>                     ACF1 Theil's U
#> Training set -0.01284696        NA
#> Test set     -0.13435655 0.9092797

MAPE:

Train = 5.81%
Test = 14.13% (overfitting)

Tindak lanjut:

Tuning alpha, beta agar menghasilkan model yang best fit. Boleh dituning secara manual, ataupun menggunakan teknik grid search yang disediakan oleh fungsi ts_grid() pada package TSstudio
Gunakan model/algoritma lain

5.2.3 Triple Exponential Smoothing / Holt-Winters

Holt winters digunakan untuk data time series yang mengandung trend dan seasonal. Holt winters melakukan smoothing terhadap komponen error/random (alpha), trend (beta), dan seasonal (gamma). Kita akan menggunakan data persentase kelahiran di New York untuk melakukan fitting model dan forecasting menggunakan metode holt winters.

birth_ts %>% 
  autoplot()

Cross validation

# 3 tahun (36 bulan) terakhir untuk test
birth_test <- tail(birth_ts, 3*12)

# sisanya untuk train
birth_train <- head(birth_ts, -3*12)

Fitting Model

Secara default HoltWinters akan menganggap data sebagai additive seasonality (seasonal = "additive"). Apabila data berupa multiplicative, gunakan parameter seasonal = "multiplicative" di dalam fungsi Holtwinters().

birth_hw <- HoltWinters(birth_train, seasonal = "additive")
birth_hw

#> Holt-Winters exponential smoothing with trend and additive seasonal component.
#> 
#> Call:
#> HoltWinters(x = birth_train, seasonal = "additive")
#> 
#> Smoothing parameters:
#>  alpha: 0.8238483
#>  beta : 0.01879443
#>  gamma: 1
#> 
#> Coefficients:
#>            [,1]
#> a   27.36749822
#> b    0.02551976
#> s1  -0.62729995
#> s2  -1.71011964
#> s3   1.64280798
#> s4   0.05574370
#> s5   1.08739310
#> s6   0.27050131
#> s7   1.54289802
#> s8   0.92194413
#> s9   0.29164988
#> s10  0.48480479
#> s11 -1.50609008
#> s12 -0.38049822

Forecasting

birth_forecast <- forecast(birth_hw, h = 3*12)

horizon = kita mau forecast berapa data ke depan?

Visualisasikan data actual dengan hasil forecasting holt winters

birth_ts %>% 
  autoplot(series = "Actual") +
  autolayer(birth_forecast$fitted, series = "Train") +
  autolayer(birth_forecast$mean, series = "Test")

Evaluasi Model

forecast::accuracy(birth_forecast, birth_test)

#>                       ME      RMSE       MAE        MPE     MAPE      MASE
#> Training set  0.09810363 0.7530131 0.5681970  0.3510623 2.347179 0.5687989
#> Test set     -0.52267504 1.0812812 0.8125558 -2.0048640 3.015850 0.8134166
#>                   ACF1 Theil's U
#> Training set 0.1490922        NA
#> Test set     0.6856166 0.6939018

MAPE:

Train: 2.347%
Test: 3.016%

5.2.4 Fungsi ets()

ets() dari library forecast yang merupakan error, trend, dan seasonal. Parameter yang digunakan, yaitu:

y: objek time series.
model: model/pola time series untuk error, trend, dan seasonal.
1. A: additive
2. M: multiplicative
3. N: none
4. Z: auto -> berdasarkan AIC (information loss, semakin kecil semakin baik)
alpha: bobot untuk smoothing error, nilainya antara 0 - 1.
beta: bobot untuk smoothing trend, nilainya antara 0 - 1.
gamma: bobot untuk smoothing seasonal, nilainya antara 0 - 1.

Jika bobot semakin mendekati 1, maka nilai prediksi/forecast lebih mempertimbangkan data terbaru.

Single Exponential Smooting/SES: ets(y, model = "*NN")

Double Exponential Smooting/Holt: ets(y, model = "**N")

Triple Exponential Smooting/Holt Winters: ets(y, model = "***")

artinya dapat diisi dengan A/M/Z

birth_ts %>% 
  autoplot()

Trend: cenderung linear -> additive Seasonality: varians konstan -> additive

birth_ets <- ets(birth_train, model = "AAA")
birth_ets

#> ETS(A,Ad,A) 
#> 
#> Call:
#>  ets(y = birth_train, model = "AAA") 
#> 
#>   Smoothing parameters:
#>     alpha = 0.7193 
#>     beta  = 0.0001 
#>     gamma = 0.0001 
#>     phi   = 0.863 
#> 
#>   Initial states:
#>     l = 26.3085 
#>     b = -0.6758 
#>     s = -0.5129 -1.2048 0.7337 0.6324 1.1251 1.526
#>            -0.0403 0.3452 -0.7701 0.9206 -2.1293 -0.6257
#> 
#>   sigma:  0.6392
#> 
#>      AIC     AICc      BIC 
#> 544.1604 550.2135 596.0508

Cobalah dengan parameter model:

ZZZ -> AIC = 544.7426 (information loss)
AAA -> AIC = 544.1604

Forecasting

birth_ets_forecast <- forecast(birth_ets, h = 3*12) # panjang data test

Visualisasikan data actual dengan hasil forecasting holt winters

birth_ts %>% 
  autoplot(series = "Actual") +
  autolayer(birth_ets_forecast$fitted, series = "Train") +
  autolayer(birth_ets_forecast$mean, series = "Test")

Evaluasi Model

forecast::accuracy(birth_ets_forecast, birth_test)

#>                       ME      RMSE       MAE        MPE     MAPE      MASE
#> Training set  0.05781606 0.5965746 0.4476847  0.1980528 1.827651 0.4481590
#> Test set     -0.05731999 0.7280613 0.5871715 -0.2885430 2.142523 0.5877936
#>                   ACF1 Theil's U
#> Training set 0.1066730        NA
#> Test set     0.5234464 0.4572179

Bandingkan performa model forecasting data birth dengan fungsi HoltWinters() dengan ets():

Fungsi HoltWinters:

Training -> MAPE 2.347%
Testing -> MAPE 3.016%

Fungsi ets (AAA):

Training -> MAPE 1.827651%
Testing -> MAPE 2.142523%

5.3 📝 SUMMARY DAY 2

TS Workflow:

Simple Moving Average (SMA):
- no trend no seasonal
- bobot untuk tiap observasi sama
Exponential smoothing:
- Single / Simple: no trend no seasonal
- Double / Holt: yes trend no seasonal
- Triple / Holt-Winters: yes trend yes seasonal
Parameter smoothing:
- alpha: untuk smoothing error
- beta: untuk smoothing trend
- gamma: untuk smoothing seasonal
- parameter smoothing semakin mendekati nilai 1 artinya bobot data baru dianggap jauh lebih penting dibandingkan data yang sudah lama
Aturan cross validation dalam TS:
- Metode pengambilan: secara berurutan
- Pengambilan data test sebaiknya data dengan waktu terbaru

DAY 3

5.4 ARIMA

Auto Regressive Integrated Moving Average (ARIMA) adalah gabungan dari dua metode, yaitu Autoregressive (AR) dan Moving Average (MA). I-nya menjelaskan Integrated. Konsep utama dari ARIMA, yaitu menjelaskan autokorelasi pada data.

5.4.1 Auto Regressive

AR(p): Autoregressive

Pada model regresi linier, target ($y$) diprediksi berdasarkan prediktor ($x$).

\[\hat{y}=\beta_0+\beta_1*x_1+...+\beta_n*x_n\]

Pada model autoregressive, target di masa depan ($\hat{y_t}$) diforecast berdasarkan nilai target di masa lalu ($y_{t-1}, ..., y_{t-n}$).

\[\hat{y}_t = c + \beta_1.y_{t-1} + \cdots + \beta_p.y_{t-p} + e_{t}\]

Dikenal istilah lag/delay, yaitu data pada waktu sebelumnya, di mana:

lag 1 = observasi pertama sebelumnya
lag 2 = observasi kedua sebelumnya
lag 3 = observasi ketiga sebelumnya, dan seterusnya

Contoh pada data harian:

lag 0 ($y_{t}$) = Kamis
lag 1 ($y_{t-1}$) = Rabu
lag 2 ($y_{t-2}$) = Selasa
lag 3 ($y_{t-3}$) = Senin

Order $p$: berapa banyak lag untuk nilai $y$ yang kita gunakan dalam melakukan forecasting

5.4.2 Moving Average

MA(q): Moving average

Nilai $\hat{y_t}$ dihasilkan dari smoothing error/random. Hasil smoothing yang dilakukan pada error/random bukan digunakan untuk melakukan forecast nilai $\hat{y_t}$, melainkan digunakan untuk mengestimasi trend.

\[\hat{y}_t = c + \theta_1.\varepsilon_{t-1} + \cdots + \theta_q.\varepsilon_{t-q} + e_{t}\]

Order $q$: berapa banyak lag untuk nilai error $\varepsilon$ yang kita gunakan dalam melakukan forecasting

ARIMA(p,d,q): non-seasonal ARIMA

\[\hat{y}_t = c + \beta_1.y_{t-1} + \cdots + \beta_p.y_{t-p} + \theta_1.\varepsilon_{t-1} + \cdots + \theta_q.\varepsilon_{t-q} + e_{t}\]

5.4.3 Integrated (Differencing)

Syarat ARIMA: time series harus stasioner (rata-rata dan variansinya konstan), artinya tidak ada trend dan tidak ada seasonal. Jika data time series belum stationer, maka perlu dilakukan differencing/Integrated pada ARIMA.

Cara yang paling umum untuk membuat data time series menjadi stasioner adalah dengan melakukan differencing diff(), yaitu mengurangi data saat ini dengan data sebelumnya (lag = 1). Terkadang, tergantung pada kompleksitas data, jumlah differencing bisa lebih dari 1 kali.

x <- c(5, 9, 3, 4, 7)
x %>% 
  diff()

#> [1]  4 -6  1  3

Pertanyaannya, mengapa kita perlu melakukan differencing untuk mendapatkan time series yang stasioner? Berikut adalah penjelasannya:

Data AirPassengers sebelum differencing memiliki korelasi sangat kuat antar lag nya:

data.frame(
  xt = c(AirPassengers),
  lag1 = lag(x = as.vector(AirPassengers), n = 1),
  lag2 = lag(x = as.vector(AirPassengers), n = 2),
  lag3 = lag(x = as.vector(AirPassengers), n = 3),
  lag4 = lag(x = as.vector(AirPassengers), n = 4),
  lag5 = lag(x = as.vector(AirPassengers), n = 5) 
) %>% 
  na.omit() %>% 
  GGally::ggcorr(label = T)

Data AirPassengers setelah differencing, korelasi antar lag tidak terlalu kuat.

data.frame(
  xt = c(diff(AirPassengers)),
  lag1 = lag(x = as.vector(diff(AirPassengers)), n = 1),
  lag2 = lag(x = as.vector(diff(AirPassengers)), n = 2),
  lag3 = lag(x = as.vector(diff(AirPassengers)), n = 3),
  lag4 = lag(x = as.vector(diff(AirPassengers)), n = 4),
  lag5 = lag(x = as.vector(diff(AirPassengers)), n = 5) 
) %>% 
  na.omit() %>% 
  GGally::ggcorr(label = T)

Ide utama melakukan differencing adalah ketika melakukan prediksi, agar multicolinearity antar lag-nya tidak terlalu besar.

5.4.4 Workflow

Import data to R

no <- read.csv("data_input/environment_1970f.csv")

Data terdiri dari 43 observasi yang mewakili kontribusi gas emisi per tahun terhadap atmosfer Indonesia (43 tahun, 1970-2012). Sebelumnya kita telah menggunakan kolom emisi CO2, sekarang mari kita gunakan informasi emisi gas Nitro Oksida (ribu metrik ton setara CO2) yaitu emisi dari pembakaran biomassa pertanian, kegiatan industri, dan pengelolaan ternak.

Exploratory data

p <- no %>% 
  ggplot(aes(x = year,
             y = Nitrous.oxide.emissions..thousand.metric.tons.of.CO2.equivalent.)) +
  geom_point() +
  geom_line() +
  labs(y = "Nitrous Oxide")

plotly::ggplotly(p) # interactive plotting

Berdasarkan line plot, apakah dapat dikatakan data sudah stasioner?

Trend: ada tapi naik tipis
Seasonality: sepertinya tidak ada pola berulang

Dari visualisasi kita belum yakin secara objektif apakah data kita memiliki trend/seasonal.

Membuat objek time series

no_ts <- ts(data = no$Nitrous.oxide.emissions..thousand.metric.tons.of.CO2.equivalent.,
            start = 1970,
            frequency = 1)

Cross-validation

# 8 tahun terakhir untuk data test, sisanya jadikan data train
no_test <- tail(no_ts, 8)
no_train <- head(no_ts, -8)

Stationarity test

ADF test

Kalau hanya dari visualisasi, penentuan apakah time series sudah stasioner akan sangat subjektif. Untuk itu, kita dapat melakukan Augmented Dickey-Fuller (ADF) test untuk menguji stationarity dengan fungsi adf.test() dari library tseries.

Hipotesis:

$H_0$: Punya unit root (tidak stasioner)
$H_1$: Tidak punya unit root (stasioner)

Harapannya adalah p-value < alpha (0.05) -> stasioner

Note: lakukan tes hanya pada data train saja

library(tseries)
adf.test(no_train)

#> 
#>  Augmented Dickey-Fuller Test
#> 
#> data:  no_train
#> Dickey-Fuller = -3.2514, Lag order = 3, p-value = 0.09523
#> alternative hypothesis: stationary

Berdasarkan ADF test, data no_train sudah stasioner atau belum?

p-value = 0.09523 alpha = 0.05

p-value > alpha, kesimpulannya gagal tolak H0 sehingga tidak stasioner

Apabila data belum stasioner, maka lakukan differencing hingga stasioner.

no_train_diff <- no_train %>% diff()

Uji lagi ADF test: harapan p-value < alpha (0.05)

no_train_diff %>% 
  adf.test()

#> 
#>  Augmented Dickey-Fuller Test
#> 
#> data:  .
#> Dickey-Fuller = -5.943, Lag order = 3, p-value = 0.01
#> alternative hypothesis: stationary

p-value = 0.01 < alpha (0.05) sehingga tolak H0

Kesimpulan: Data no_train perlu dilakukan differencing sebanyak satu kali hingga stasioner.

Mari kita bandingkan no_train tanpa dan dengan differencing melalui visualisasi:

no_train %>% 
  autoplot(series = "Actual") +
  autolayer(no_train_diff, series = "Diff 1")

ARIMA(p,d,q) = AR(p) I(d) MA(q)

Banyaknya differencing yang dilakukan sebelumnya menjadi bagian order d untuk komponen Integrated I(d), sehingga menjadi model ARIMA(p,1,q)

Fitting model otomatis menggunakan auto.arima()

PENTING: data yang di-fit pada fungsi auto.arima() adalah data sebelum differencing

Gunakan parameter seasonal = F untuk ARIMA biasa

no_auto <- auto.arima(y = no_train, seasonal = F)
summary(no_auto)

#> Series: no_train 
#> ARIMA(0,1,1) 
#> 
#> Coefficients:
#>           ma1
#>       -0.8150
#> s.e.   0.0919
#> 
#> sigma^2 estimated as 3317659997:  log likelihood=-420.96
#> AIC=845.93   AICc=846.32   BIC=848.98
#> 
#> Training set error measures:
#>                    ME     RMSE      MAE      MPE     MAPE      MASE        ACF1
#> Training set 10972.01 55929.24 31826.05 0.950838 23.63757 0.7418046 -0.04361855

Apa artinya ARIMA(0,1,1)?

I(d): d = 1 -> artinya data awal perlu di differencing 1x sebelum dilakukan pemodelan
AR(p): p = 0 -> artinya model tidak menggunakan data lag sama sekali
MA(q): q = 1 -> artinya model menggunakan error sampai lag 1 saja

Fitting model manual dengan Arima() dari package forecast

# membuat model dengan order p dan q nya asal
# misal ARIMA(1,1,2)
no_manual <- Arima(y = no_train, order = c(1,1,2))
summary(no_manual)

#> Series: no_train 
#> ARIMA(1,1,2) 
#> 
#> Coefficients:
#>           ar1     ma1      ma2
#>       -0.8050  0.1723  -0.8276
#> s.e.   0.1183  0.1683   0.1480
#> 
#> sigma^2 estimated as 3169988254:  log likelihood=-419.85
#> AIC=847.7   AICc=849.08   BIC=853.8
#> 
#> Training set error measures:
#>                    ME     RMSE      MAE      MPE     MAPE      MASE      ACF1
#> Training set 10336.83 52987.77 29322.74 1.500105 21.51614 0.6834571 -0.119463

Bandingkan kedua performa model:

auto: ARIMA(0,1,1) -> MAPE 23.63757%
manual: ARIMA(1,1,2) -> MAPE 21.51614%

Silahkan jadikan model auto sebuah baseline model, tapi perlu ingat bahwa model auto belum tentu menjadi model yang terbaik

5.4.5 Menentukan p dan q secara manual

Menentukan ARIMA (p, d, q) secara manual

Langkah 1: Tentukan d (dari differencing)

Langkah 2: Tentukan p dan q

tribble(~model, ~ACF, ~PACF,
        "AR(p)", "Die Down", "Cut off after lag p",
        "MA(q)", "Cut off after lag q", "Die Down",
        "ARMA(p, q)", "Die Down", "Die Down",
        "No AR(p) or MA(q)", "No spike", "No spike"
        )

Perbedaan PACF dan ACF:

ACF (Auto Correlation Function) digunakan untuk menentukan order q MA(q)
PACF (Partial Auto Correlation Function) digunakan untuk menentukan order p AR(p)

Misal hari ini Kamis:

ACF: memperhitungkan korelasi yang transitif, contoh:
- ACF lag 1: Korelasi data Kamis dengan data Rabu
- ACF lag 2: Korelasi data Kamis dengan data Selasa, tapi masih memperhitungkan korelasi Kamis-Rabu dan Rabu-Selasa
PACF: tidak memperhitungkan korelasi yang transitif, hanya korelasi parsialnya, contoh:
- PACF lag 1: Korelasi data Kamis dengan data Rabu
- PACF lag 2: Korelasi data Kamis dengan data Selasa, tanpa mempertimbangkan efek korelasi pada Rabu

Die down vs Cut off

Die down/decay: mengalami penurunan atau kenaikan secara lambat
Cut off: mengalami penurunan atau kenaikan secara cepat

Mari kita cek ACF dan PACF untuk no_train_diff

PENTING: Saat plot ACF dan PACF, gunakan data yang sudah di-differencing (karena kita akan melakukannya secara manual)

# fungsi tsdisplay dari package forecast
tsdisplay(no_train_diff)

Melakukan plot ACF dan PACF bisa menggunakan tsdisplay() namun agar tampil lebih besar, kita buat custom function plot_cf() sebagai berikut:

# opsional: custom function
library(patchwork)

plot_cf <- function(x, lag.max = NULL){
  p1 <- ggAcf(x, lag.max=lag.max) + labs(title = "Plot ACF and PACF")
  p2 <- ggPacf(x, lag.max=lag.max) + labs(title = NULL)
  p1 / p2
}

plot_cf(no_train_diff)

Garis putus-putus horizontal adalah batasan signifikansi autokorelasi, range nilainya adalah +- 1.96 / akar(T), dengan T adalah panjang time series.

Dari plot ACF-PACF dapat disimpulkan:

ACF: cut off di lag = 1
PACF: cut off di lag = 4

Setelah itu mari tentukan order:

p = 1, 2, 4 (signifikan)
d = 1
q = 1

Kombinasi model ARIMA(p,d,q):

ARIMA(1,1,1)
ARIMA(2,1,1)
ARIMA(4,1,1)

Istilah: Spike -> lag dengan nilai autokorelasi yang signifikan

Fitting kombinasi model

Note: fit data yang sebelum differencing

no_arima1 <- Arima(y = no_train, order = c(1,1,1))
no_arima2 <- Arima(y = no_train, order = c(2,1,1))
no_arima3 <- Arima(y = no_train, order = c(4,1,1))

Cek error di training set dengan summary()

summary(no_arima1)

#> Series: no_train 
#> ARIMA(1,1,1) 
#> 
#> Coefficients:
#>          ar1      ma1
#>       0.0267  -0.8216
#> s.e.  0.1943   0.1019
#> 
#> sigma^2 estimated as 3420512259:  log likelihood=-420.96
#> AIC=847.91   AICc=848.71   BIC=852.49
#> 
#> Training set error measures:
#>                    ME    RMSE      MAE     MPE     MAPE      MASE        ACF1
#> Training set 11073.59 55922.5 31760.79 1.07617 23.49991 0.7402834 -0.06073625

summary(no_arima2)

#> Series: no_train 
#> ARIMA(2,1,1) 
#> 
#> Coefficients:
#>           ar1      ar2      ma1
#>       -0.0445  -0.2210  -0.7471
#> s.e.   0.2000   0.1822   0.1416
#> 
#> sigma^2 estimated as 3379874911:  log likelihood=-420.28
#> AIC=848.55   AICc=849.93   BIC=854.66
#> 
#> Training set error measures:
#>                    ME     RMSE      MAE       MPE     MAPE      MASE
#> Training set 10008.72 54713.83 30326.12 0.6124584 22.41164 0.7068441
#>                     ACF1
#> Training set -0.06979143

summary(no_arima3)

#> Series: no_train 
#> ARIMA(4,1,1) 
#> 
#> Coefficients:
#>           ar1      ar2      ar3      ar4      ma1
#>       -0.5841  -0.6503  -0.4078  -0.3948  -0.1906
#> s.e.   0.3395   0.2617   0.2346   0.1698   0.3816
#> 
#> sigma^2 estimated as 3257345416:  log likelihood=-418.67
#> AIC=849.34   AICc=852.45   BIC=858.5
#> 
#> Training set error measures:
#>                    ME     RMSE      MAE       MPE     MAPE      MASE
#> Training set 6516.114 51951.36 30469.81 -1.359721 23.11542 0.7101932
#>                     ACF1
#> Training set -0.03541161

OPSIONAL: Custom function compare_fitting() untuk membandingkan performa di data train agar tampilan lebih rapi

# opsional: custom function
compare_fitting <- function(x){
  lapply(x, function(x) summary(x)["Training set", ]) %>% 
    lapply(., function(x) x %>% t() %>% as.data.frame) %>% 
    bind_rows() %>% 
    mutate(model = names(x)) %>%
    select(model, everything())
}

model_list <- list(
  "ARIMA 011" = no_auto,
  "ARIMA 111" = no_arima1,
  "ARIMA 211" = no_arima2,
  "ARIMA 411" = no_arima3)

compare_fitting(model_list)

#> Series: no_train 
#> ARIMA(0,1,1) 
#> 
#> Coefficients:
#>           ma1
#>       -0.8150
#> s.e.   0.0919
#> 
#> sigma^2 estimated as 3317659997:  log likelihood=-420.96
#> AIC=845.93   AICc=846.32   BIC=848.98
#> 
#> Training set error measures:
#>                    ME     RMSE      MAE      MPE     MAPE      MASE        ACF1
#> Training set 10972.01 55929.24 31826.05 0.950838 23.63757 0.7418046 -0.04361855
#> Series: no_train 
#> ARIMA(1,1,1) 
#> 
#> Coefficients:
#>          ar1      ma1
#>       0.0267  -0.8216
#> s.e.  0.1943   0.1019
#> 
#> sigma^2 estimated as 3420512259:  log likelihood=-420.96
#> AIC=847.91   AICc=848.71   BIC=852.49
#> 
#> Training set error measures:
#>                    ME    RMSE      MAE     MPE     MAPE      MASE        ACF1
#> Training set 11073.59 55922.5 31760.79 1.07617 23.49991 0.7402834 -0.06073625
#> Series: no_train 
#> ARIMA(2,1,1) 
#> 
#> Coefficients:
#>           ar1      ar2      ma1
#>       -0.0445  -0.2210  -0.7471
#> s.e.   0.2000   0.1822   0.1416
#> 
#> sigma^2 estimated as 3379874911:  log likelihood=-420.28
#> AIC=848.55   AICc=849.93   BIC=854.66
#> 
#> Training set error measures:
#>                    ME     RMSE      MAE       MPE     MAPE      MASE
#> Training set 10008.72 54713.83 30326.12 0.6124584 22.41164 0.7068441
#>                     ACF1
#> Training set -0.06979143
#> Series: no_train 
#> ARIMA(4,1,1) 
#> 
#> Coefficients:
#>           ar1      ar2      ar3      ar4      ma1
#>       -0.5841  -0.6503  -0.4078  -0.3948  -0.1906
#> s.e.   0.3395   0.2617   0.2346   0.1698   0.3816
#> 
#> sigma^2 estimated as 3257345416:  log likelihood=-418.67
#> AIC=849.34   AICc=852.45   BIC=858.5
#> 
#> Training set error measures:
#>                    ME     RMSE      MAE       MPE     MAPE      MASE
#> Training set 6516.114 51951.36 30469.81 -1.359721 23.11542 0.7101932
#>                     ACF1
#> Training set -0.03541161

Forecasting

Lakukan forecasting untuk semua kombinasi model

no_auto_f <- forecast(no_auto, h = 8) # panjang data test
no_arima1_f <- forecast(no_arima1, h = 8)
no_arima2_f <- forecast(no_arima2, h = 8)
no_arima3_f <- forecast(no_arima3, h = 8)

Evaluate model

OPSIONAL: Custom function compare_forecast() untuk membandingkan performa di data test dengan tampilan lebih rapi

# opsional: custom function
compare_forecast <- function(x, test){
  lapply(x, function(x) forecast::accuracy(x, test)["Test set", ]) %>%
    lapply(., function(x) x %>% t() %>% as.data.frame) %>% 
    bind_rows() %>% 
    mutate(model = names(x)) %>%
    select(model, everything())
}

forecast_list <- list(
  "ARIMA 011" = no_auto_f,
  "ARIMA 111" = no_arima1_f,
  "ARIMA 211" = no_arima2_f,
  "ARIMA 411" = no_arima3_f)

compare_forecast(forecast_list, no_test)

Visualisasikan hasil forecasting model

no_ts %>% 
  autoplot(series = "Actual") +
  autolayer(no_auto_f$mean, series = "ARIMA(0,1,1)") +
  autolayer(no_arima1_f$mean, series = "ARIMA(1,1,1)") +
  autolayer(no_arima2_f$mean, series = "ARIMA(2,1,1)") +
  autolayer(no_arima3_f$mean, series = "ARIMA(4,1,1)")

Kesimpulan:

Model terbaik di data test adalah ARIMA(4,1,1)
Training MAPE = 23.11542% sedangkan Testing MAPE = 28.16403%
Apakah underfit/best fit/overfit? underfitting (dengan asumsi kebutuhan kita di bawah 10%)

5.5 📝 SUMMARY DAY 3

ARIMA(p,d,q)
- Auto Regressive AR(p): regresi untuk data lag nya, sebanyak p
- Integrated I(d): differencing sebanyak d kali agar datanya stasioner
- Moving Average MA(q): regresi untuk error lag nya, sebanyak q
ARIMA Workflow:
1. Data preprocessing time series, memastikan sudah dalam bentuk time series object ts()
2. Cross-validation: train untuk modeling, test untuk evaluation
3. Uji stasioner: dengan adf.test(), harapannya p-value < alpha agar data stasioner
4. Apabila belum stasioner, maka lakukan differencing hingga stasioner
5. Buat baseline model dengan auto.arima(obj_ts, seasonal=F)
6. Tuning model dengan menentukan order p dan q secara manual
- Order p ditentukan dari plot PACF
- Order q ditentukan dari plot ACF
- Memperhatikan pola die down dan cut-off dari plot tersebut
1. Fitting model dari kombinasi order p,d,q nya
2. Forecasting sesuai dengan panjang data test
3. Evaluasi model:
- Membandingkan error di data test
- Visualisasi hasil forecasting

DAY 4

5.6 SARIMA (Seasonal ARIMA)

Seasonal ARIMA adalah metode ARIMA di mana data time series memiliki pola seasonal. Tahapan dalam melakukan pemodelan menggunakan SARIMA sama seperti saat membuat pemodelan ARIMA.

ARIMA: ARIMA(p,d,q)
SARIMA: ARIMA(p,d,q)(P,D,Q)[m]

Keterangan:

(p,d,q) sebagai index utama (cara penentuan sama seperti ARIMA biasa)
(P,D,Q) sebagai index seasonal -> lag pada kelipatan nilai m (frequency)
[m] merupakan nilai frequency pada time series

Jika kita menggunakan model SARIMA, diperlukan langkah tambahan ketika melakukan differencing, yaitu differencing sesuai pola berulang/frequency untuk menghilangkan seasonality.

5.6.1 Workflow

Data DAUTONSA.csv merupakan penjualan motor di suatu dealer dari tahun 1967-2019.

DATE: tanggal saat dilakukan pencatatan penjualan
DAUTONSA: rata-rata penjualan bulanan

sales <- read.csv("data_input/DAUTONSA.csv")
glimpse(sales)

#> Rows: 102
#> Columns: 2
#> $ DATE     <chr> "2011-01-01", "2011-02-01", "2011-03-01", "2011-04-01", "2011~
#> $ DAUTONSA <dbl> 258.576, 338.615, 445.315, 405.418, 362.214, 348.950, 329.517~

Cek apakah data sudah memenuhi time series yang baik?

sales %>% 
  mutate(DATE = ymd(DATE)) %>% # ubah jadi tipe data date
  arrange(DATE) %>% # mengurutkan berdasarkan tanggal
  pad() %>% # menambal tanggal (padding)
  anyNA()

#> [1] FALSE

Artinya data kita sudah memenuhi karakteristik data time series yang baik

Data wrangling

Ubah DATE menjadi tipe data tanggal
Ekstrak informasi bulan pada DATE untuk kebutuhan visualisasi

sales <- sales %>% 
  mutate(DATE = ymd(DATE),
         MONTH = month(DATE, label = T))

Visualisasi

sales %>% 
  ggplot(aes(DATE, DAUTONSA)) +
  geom_point() +
  geom_point(sales %>% filter(MONTH == "Jan"), 
             mapping = aes(DATE, DAUTONSA), col = "blue") +
  geom_line()

Membuat object time series

sales_ts <- ts(data = sales$DAUTONSA,
               start = 2011,
               frequency = 12) # data bulanan, pola tahunan

Decompose

sales_ts %>% 
  decompose() %>% 
  autoplot()

Cross-validation

# 1.5 tahun untuk test
sales_test <- tail(sales_ts, 1.5*12)

# 7 tahun untuk train, sisanya
sales_train <- head(sales_ts, -1.5*12)

Cek stationarity data

Gunakan adf.test(): - H0 = Tidak Stationary - H1 = Stationary

sales_train %>% 
  adf.test()

#> 
#>  Augmented Dickey-Fuller Test
#> 
#> data:  .
#> Dickey-Fuller = -2.9359, Lag order = 4, p-value = 0.192
#> alternative hypothesis: stationary

Kesimpulan: p-value > alpha (gagal tolak H0), artinya data sales_train masih tidak stasioner

ARIMA: diff(lag=1)

Feb 2021 dikurangkan dengan Jan 2021
Jan 2021 dikurangkan dengan Des 2020
dst

SARIMA: diff(lag=frequency)

Kasus frequency = 12:

Feb 2021 dikurangkan dengan Feb 2020
Jan 2021 dikurangkan dengan Jan 2020
dst

sales_train_diff <- sales_train %>% 
  diff(lag=12) %>% # differencing untuk komponen seasonal
  diff(lag=1) # differencing untuk trend

sales_train_diff %>% 
  adf.test()

#> 
#>  Augmented Dickey-Fuller Test
#> 
#> data:  .
#> Dickey-Fuller = -4.7557, Lag order = 4, p-value = 0.01
#> alternative hypothesis: stationary

Kesimpulan:

d = 1 -> melakukan differencing 1x untuk komponen trend
D = 1 -> melakukan differencing 1x untuk komponen seasonalnya

Note: lakukan differencing seasonal terlebih dahulu, baru trend-nya agar komponen waktu tidak bergeser. Urutannya:

Cek dahulu apakah data stasioner/tidak? dengan ADF test
Apabila tidak stasioner, lakukan differencing SEASONAL terlebih dahulu lalu cek ADF test lagi
Apabila masih tidak stasioner, lakukan differencing TREND (lag = 1) lalu cek ADF test lagi
Apabila masih tidak stasioner, lakukan differencing SEASONAL (tapi letakkan setelah SEASONAL yang di step 2) -> lalu cek ADF test lagi
Apabila masih tidak stasioner, lakukan differencing TREND (diletakkan setelah TREND yang di step 3) lalu cek ADF test lagi
Ulangi langkah sebelumnya hingga stasioner

Baseline model

Gunakan auto.arima(), namun set parameter seasonal = T agar menjadi auto SARIMA

ARIMA: seasonal = F SARIMA: seasonal = T

sales_auto <- auto.arima(sales_train, seasonal = T)
sales_auto

#> Series: sales_train 
#> ARIMA(0,1,1)(0,1,0)[12] 
#> 
#> Coefficients:
#>           ma1
#>       -0.6468
#> s.e.   0.0900
#> 
#> sigma^2 estimated as 892.7:  log likelihood=-341.71
#> AIC=687.41   AICc=687.59   BIC=691.94

Auto SARIMA: ARIMA(0,1,1)(0,1,0)[12]

Visualisasi hasil fitting model

sales_train %>%
  autoplot(series = "Actual") +
  autolayer(sales_auto$fitted, series = "ARIMA(0,1,1)(0,1,0)[12]")

Manual SARIMA

Gunakan custom function plot_cf() untuk menentukan order p, q, P, Q

# visualisasi PACF ACF
plot_cf(sales_train_diff, lag.max = 50)

PACF

non-seasonal: p = 1,2
seasonal (freq = 12): P = 0

ACF

non-seasonal: q = 1
seasonal (freq = 12): Q = 0

Kombinasi model ARIMA(p,d,q)(P,D,Q)[m]:

ARIMA(1,1,1)(0,1,0)[12]
ARIMA(2,1,1)(0,1,0)[12]

Fitting model manual

Parameter:

non-seasonal: order = c(p,d,q)
seasonal: seasonal = c(P,D,Q)

sales_sarima1 <- Arima(sales_train,
                       order = c(1,1,1),
                       seasonal = c(0,1,0))

Gunakan custom function compare_fitting() untuk membandingkan performa model di data train

model_list <- list(
  "ARIMA(0,1,1)(0,1,0)[12]" = sales_auto,
  "ARIMA(1,1,1)(0,1,0)[12]" = sales_sarima1)

compare_fitting(model_list)

#> Series: sales_train 
#> ARIMA(0,1,1)(0,1,0)[12] 
#> 
#> Coefficients:
#>           ma1
#>       -0.6468
#> s.e.   0.0900
#> 
#> sigma^2 estimated as 892.7:  log likelihood=-341.71
#> AIC=687.41   AICc=687.59   BIC=691.94
#> 
#> Training set error measures:
#>                    ME     RMSE      MAE      MPE     MAPE      MASE        ACF1
#> Training set -4.11363 27.27477 19.30697 -1.13262 4.466261 0.4645816 -0.02948427
#> Series: sales_train 
#> ARIMA(1,1,1)(0,1,0)[12] 
#> 
#> Coefficients:
#>           ar1      ma1
#>       -0.0095  -0.6415
#> s.e.   0.1808   0.1361
#> 
#> sigma^2 estimated as 905.6:  log likelihood=-341.71
#> AIC=689.41   AICc=689.77   BIC=696.2
#> 
#> Training set error measures:
#>                     ME     RMSE      MAE       MPE     MAPE      MASE
#> Training set -4.092282 27.27413 19.30068 -1.128451 4.464672 0.4644302
#>                     ACF1
#> Training set -0.02586758

Forecasting

sales_auto_f <- forecast(sales_auto, h = 1.5*12)
sales_sarima1_f <- forecast(sales_sarima1, h = 1.5*12)

Evaluasi model

Gunakan custom function compare_forecast() untuk membandingkan performa model di data test

forecast_list <- list(
  "ARIMA(0,1,1)(0,1,0)[12]" = sales_auto_f,
  "ARIMA(1,1,1)(0,1,0)[12]" = sales_sarima1_f)

compare_forecast(forecast_list, sales_test)

Visualisasi hasil forecasting

sales_ts %>% 
  autoplot(series = "Actual") +
  autolayer(sales_auto_f$mean, series = "ARIMA(0,1,1)(0,1,0)[12]") +
  autolayer(sales_sarima1_f$mean, series = "ARIMA(1,1,1)(0,1,0)[12]")

Optional: untuk visualisasi hasil forecast yang interaktif menggunakan fungsi test_forecast() dari package TSstudio

test_forecast(actual = sales_ts, 
              forecast.obj = sales_auto_f, 
              train = sales_train, 
              test = sales_test)

6 Assumption

Asumsi pada time series diujikan untuk mengukur apakah residual yang peroleh dari hasil modeling sudah cukup baik untuk menggambarkan dan menangkap informasi pada data. Mengapa menggunakan residual data? Karena kita dapat mendapatkan informasi dari data aktual maupun dari hasil prediksi menggunakan model.

Metode forecasting yang baik menghasilkan nilai residual berikut ini:

Residual yang tidak berkorelasi. Apabila terdapat residual yang berkorelasi, artinya masih terdapat informasi yang tertinggal yang seharusnya digunakan untuk menghitung hasil forecast.
Residual berdistribusi normal di sekitaran nilai 0.

Untuk memastikan bahwa residual yang dihasilkan memiliki kriteria diatas, maka terdapat asumsi untuk mengujinya.

Note: Cek asumsi cukup dilakukan pada model akhir yang akan dipilih

6.1 No auto-correlation residual

Untuk cek ada/tidaknya autokorelasi pada hasil forecasting dapat menggunakan visualisasi plot residual ACF dengan fungsi Acf(model$residuals)

# $residuals selisih nilai actual dan forecast di data train
Acf(sales_sarima1$residuals)

Berdasarkan plot ACF di atas, pada lag 8 dan 14 terdapat autokorelasi. Namun apakah secara keseluruhan lag, autokorelasi tersebut signifikan? Visualisasi bersifat subjektif, apabila ingin lebih pasti, lakukan uji Ljung-box menggunakan fungsi Box.test(model$residuals, type = "Ljung-Box") dengan hipotesis:

$H_0$: tidak terdapat autokorelasi pada residual (autokorelasi = 0)
$H_1$: terdapat autokorelasi pada residual (autokorelasi != 0)

Harapan: p-value > alpha (0.05)

Box.test(sales_sarima1$residuals, type = "Ljung-Box")

#> 
#>  Box-Ljung test
#> 
#> data:  sales_sarima1$residuals
#> X-squared = 0.058239, df = 1, p-value = 0.8093

Kesimpulan: p-value > alpha maka gagal tolak H0, sehingga tidak terdapat autokorelasi pada residual

6.2 Normality residual

Uji ini sama persis dengan uji normalitas pada course Regression Model.

hist(sales_sarima1$residuals, breaks = 15)

Berdasarkan histogram di atas, distribusi dari residual masih rancu apakah normal atau tidak. Apabila ingin lebih pasti, lakukan uji shapiro test menggunakan fungsi shapiro.test(model$residuals) dengan hipotesis:

$H_0$: residual menyebar normal
$H_1$: residual tidak menyebar normal

Harapan: p-value > alpha (0.05)

shapiro.test(sales_sarima1$residuals)

#> 
#>  Shapiro-Wilk normality test
#> 
#> data:  sales_sarima1$residuals
#> W = 0.96604, p-value = 0.02575

Kesimpulan: p-value < alpha sehingga tolak H0, maka residual tidak menyebar normal (tidak memenuhi asumsi)

Keep in mind: model yang kita bangun masih bias

7 (Additional) STLM

STLM = Seasonal Trend with Loess Model

Apabila dalam decompose biasa, dalam mendapatkan komponen trend dengan cara central moving average (CMA) di mana secara konsep setiap data yang ingin dirata-ratakan diberikan bobot yang sama sesuai ordo yang ditetapkan. Karena merata-ratakan data tengahnya hasilnya kita kehilangan data awal dan data akhirnya, sehingga ada beberapa informasi yang hilang.

Ada salah satu cara untuk mendapatkan decompose data namun tetap mempertahankan informasi dari seluruh data yang kita miliki yaitu dengan menggunakan STL (Seasonal Trend with Loess). STL secara konsep akan melakukan smoothing terhadap data tetangga setiap masing-masing observasi dengan memberikan bobot yang lebih berat terhadap data yang dekat dengan observed data.

Kekurangan dari STL hanya bisa melakukan decompose pada additive data, apabila terdapat multiplicative data dapat menggunakan transformasi log().

Untuk memodelkan hasil STL, kita juga dapat menerapkan metode exponential smoothing (ETS) dan ARIMA. Selain itu, STLM dapat digunakan sebagai alternative cara untuk menangkap seasonal yang belum bisa ditangkap oleh metode ETS dan ARIMA biasa.

Parameter stlm():

y : object time series
s.window : pola seasonal yang ingin ditangkap (frequency)
method : metode forecast yang akan digunakan, tersedia ets dan arima

sales_stlm <- stlm(y = sales_train,
                   s.window = 12,
                   method = "arima")

summary(sales_stlm$model)

#> Series: x 
#> ARIMA(2,1,1) 
#> 
#> Coefficients:
#>           ar1      ar2    ma1
#>       -1.3573  -0.5765  0.970
#> s.e.   0.0903   0.0897  0.046
#> 
#> sigma^2 estimated as 562.9:  log likelihood=-380.07
#> AIC=768.15   AICc=768.66   BIC=777.82
#> 
#> Training set error measures:
#>                     ME     RMSE      MAE         MPE     MAPE      MASE
#> Training set 0.4646435 23.15404 18.12287 -0.08624959 4.369887 0.4372299
#>                     ACF1
#> Training set -0.01882036

Berikut adalah hasil decompose dibalik STLM:

Perhatikan panjang trend sama dengan panjang data train (tidak ada yang NA seperti decompose tradisional)
Perhatikan seasonal tiap periodenya dapat berubah

# decompose biasa
sales_train %>% 
  decompose() %>% 
  autoplot()

# decompose stl
sales_stlm$stl %>% 
  autoplot()

Forecasting

sales_stlm_forecast <- forecast(sales_stlm, h = 18)
forecast::accuracy(sales_stlm_forecast$mean, sales_test)

#>                 ME     RMSE      MAE       MPE     MAPE      ACF1 Theil's U
#> Test set -37.34721 46.99031 38.06995 -11.28805 11.52644 0.4927543 0.9266351

Visualisasi

sales_ts %>% 
  autoplot(series = "Actual") +
  autolayer(sales_stlm_forecast$fitted, series = "Train") +
  autolayer(sales_stlm_forecast$mean, series = "Test")

8 📝 SUMMARY DAY 4

Seasonal ARIMA: ARIMA(p,d,q)(P,D,Q)[m]
- p,d,q sama persis dengan ARIMA biasa
- m adalah frequency seasonal time series
- D adalah banyaknya differencing pada lag m
- P dan Q ditentukan dari PACF dan ACF tapi hanya pada kelipatan lag m saja
Workflow SARIMA mirip dengan workflow ARIMA tetapi berbeda di tahap berikut:
- Lakukan differencing terhadap seasonalnya terlebih dahulu untuk mendapatkan order D, kemudian terhadap trend untuk order d
- Buat baseline model dengan auto.arima(obj_ts, seasonal=T)
- Selain menentukan p dan q, tentukan juga P dan Q
Asumsi model Time Series Forecasting:
- No auto-correlation residual
- Normality residual

Time Series: In-class Materials

Tomy Tjandra

May 28, 2021

1 Pendahuluan

1.1 Mind Map

1.2 Library

2 Time Series and Forecasting

3 Time Series Data

3.1 Characteristics

3.2 Time Series Object

3.3 Exploring nybirth Data

3.4 Decomposition

3.5 Additive and Multiplicative

3.5.1 (Optional) Perhitungan manual decompose Additive

3.5.2 Multiplicative Time Series

3.5.3 (Opsional) Perhitungan manual decompose Multiplicative

3.6 Seasonality Analysis

3.7 Seasonality Adjusted

4 📝 SUMMARY DAY 1

5 Forecasting

5.1 Simple Moving Average (SMA)

5.2 Exponential Smoothing

5.2.1 Simple Exponential Smoothing (SES)

5.2.2 Double Exponential Smoothing / Holt Exponential Smoothing

5.2.3 Triple Exponential Smoothing / Holt-Winters

5.2.4 Fungsi ets()

5.3 📝 SUMMARY DAY 2

5.4 ARIMA

5.4.1 Auto Regressive

5.4.2 Moving Average

5.4.3 Integrated (Differencing)

5.4.4 Workflow

5.4.5 Menentukan p dan q secara manual

5.5 📝 SUMMARY DAY 3

5.6 SARIMA (Seasonal ARIMA)

5.6.1 Workflow

6 Assumption

6.1 No auto-correlation residual

6.2 Normality residual

7 (Additional) STLM

8 📝 SUMMARY DAY 4

9 References

9.1 Supporting Material

9.2 TS Forecasting Model

3.3 Exploring `nybirth` Data